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ABSTRACT 

Orbital motion in triaxial nuclei with central point masses, representing 
supermassive black holes, is investigated. The stellar density is assumed to 
follow a power law, p oc r~ 7 , with 7 = 1 or 7 = 2. At low energies the motion is 
essentially regular; the major families of orbits are the tubes and the pyramids. 
Pyramid orbits are similar to box orbits but have their major elongation parallel 
to the short axis of the figure. A number of regular orbit families associated 
with resonances also exist, most prominently the banana orbits, which are also 
elongated parallel to the short axis. At a radius where the enclosed stellar mass 
is a few times the black hole mass, the pyramid orbits become stochastic. The 
energy of transition to this "zone of chaos" is computed as a function of 7 and 
of the shape of the stellar figure; it occurs at lower energies in more elongated 
potentials. Our results suggest that supermassive black holes may place tight 
constraints on departures from axisymmetry in galactic nuclei, both by limiting 
the allowed shapes of regular orbits and by inducing chaos. 
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1. Introduction 

An important change in our thinking about galaxy dynamics took place during 
the last decade, when it was recognized that the central densities of early-type galaxies 
and spheroids are generically very high. Evidence for large central masses came from 
high-resolution kinematical studies of nuclear stars and gas, which revealed the presence 
in roughly a dozen galaxies of compact dark objects with masses 1O 6 ' 5_9 ' 5 M , presumably 
supermassive black holes ( [Ford et al. 1998 ). Observations with HST also demonstrated very 



high stellar densities at the centers of early-type galaxies QCrane et al. 199~3| ; [Ferrarese et 



al. 1994|). Low luminosity ellipticals have density profiles that increase as unbroken power 



laws at small radii, p ~ r -7 , with 7^2. But Kormendy pointed out already in 1985 that 
bright galaxies also have central brightness profiles that deviate systematically from that of 
an isothermal core. The significance of this deviation was not recognized for ten more years 
due to an optical illusion associated with projection onto the plane of the sky. A luminosity 
density that varies as r -7 at small radii generates a power-law cusp in projection only if 
7 > 1. When 7 = 1, the central surface brightness is logarithmically divergent (e.g. Pehnen| 



1993| , Fig. 1), and the surface brightness profile differs only subtly from that of a galaxy 
with an isothermal core. Nonparametric deprojection of the luminosity profiles of bright 
galaxies flMerritt fc Fridman 1995| ; |Gebhardt et al. 1996|) revealed that they too harbor 
power-law cusps but with indices 7 < 1. 

While the origin of the power-law cusps is not clear, there are hints that they may be 
associated with black holes. Steep cusps, 7 ~ 2, form naturally in stellar systems where the 
black holes grow on time scales long compared to crossing times (Peebles 1972 ; Quinlan, 



[Hcrnquist fc Sigurdsson 1995| ). No universally accepted model has yet been proposed for 
the origin of the weak cusps, but we note here one feature that suggests a link to black 
holes. Two galaxies, NGC 3379 and M87, have weak cusps with well-determined structural 
parameters and also have black holes with accurately determined masses. Table 1 gives for 
each galaxy the "break" radius r& at which the central power law turns over to a steeper 
outer profile; and also the stellar mass M* contained within r^. In both galaxies, M*(r&) is 
identical within the uncertainties with M., the mass of the black hole; this is in spite of 
a factor of ~ 30 difference in M, and ~ 6 in r&. This rough equality, and the exclusive 
association of weak cusps with bright galaxies, is consistent with models in which weak 
cusps are produced following galaxy mergers by the ejection of stars from the nucleus 
by binary black holes ( Ebisuzaki, Makino fc Okumura 1991| ; |Makino 1997| ; |Quinlan~%; 
[Hernquist 19971) . 



In a fixed stellar potential, the gravitational influence of the black hole is limited to 
stars with pericenters r p < r g where M*(r g ) = M,. In an axisymmetric galaxy, r p is bounded 
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from below by the orbital angular momentum L z about the symmetry axis and only a small 
fraction of the stars, most of which are confined to the nucleus, can be strongly affected 
by the black hole. But the gravitational influence of a black hole can extend far beyond 
the nucleus in a non-axisymmetric galaxy since orbital angular momenta are not conserved 
and stars with arbitrarily large energies can pass close to the center flGerhard fc Binney 



|1985|) . In a triaxial potential containing a central point mass, the phase space divides 
naturally into three regions depending on distance from the center. In the innermost region, 
t ~ T gi the potential is dominated by the black hole and the motion is essentially regular. 
Sridhar & Touma (1999) and Sambhus & Sridhar (2000) demonstrated the regularity of the 
motion in black hole nuclei with constant stellar densities, and Merritt & Valluri (1999) 
found similar results for motion in triaxial models with weak cusps, 7 = 0.5. Here we 
extend those results to cusps with the steeper power laws characteristic of most galaxies. 
At intermediate radii, the black hole acts as a scattering center rendering almost all of the 
center-filling orbits stochastic. This "zone of chaos" extends outward from a few times r g 
to a radius where the enclosed stellar mass is roughly 10 2 times the mass of the black hole 
( [Merritt 1999|) . In the outermost region, the phase space is a complex mixture of chaotic 



and regular trajectories, including resonant box orbits that remain stable by avoiding the 
center ( Carpintero fc Aguilar 19"98| ; Papaphilippou fc Laskar 1998 ; Valluri fc Merritt 1998 ; 



Wachlin & Ferraz-Mello 1998 



The focus of the present study is on the inner two regions and specifically on the 
transition from ordered motion near the black hole to chaotic motion at r > r g . As noted 
above, the break radius is of order r g in the two bright elliptical galaxies where both 
radii can be accurately measured. The sudden change in the orbital behavior near r g might 
therefore imply a change in the three-dimensional shapes of galaxies near rv In fact there 
is some evidence for changes in ellipticity and boxiness in bright elliptical galaxies at r « rj 
( [Ryden 1999| ; |Quillen 1999| ; [Bender fc Saglia 1999| ). The work presented here is a prelude 



to full self-consistency studies which will place more rigorous constraints on the allowed 
shapes of triaxial black-hole nuclei. 

In §2 we present our model for the stellar density and calculate the gravitational 
potential and forces. The families of orbits are discussed in §3, and the transition from 
regular to chaotic motion is discussed in §4. §5 sums up. 
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2. Mass Model, Potential, Forces 

We model the stellar nucleus as a triaxial spheroid with a power-law dependence of 
density on radius: 

p* = p m" 7 , (la) 
a" 1 er c z 

within a bounding ellipsoid m = m max . The isodensity surfaces are concentric ellipsoids 
with fixed axis ratios a : b : c; without loss of generality, we assume a > b > c. Triaxiality is 
measured by the quantity T which is defined as 

a 2 - b 2 

T = (2) 
er — er 

so that prolate galaxies (b = c) have T = 1 and oblate galaxies (a = b) have T = 0. Since 
the models are scale-free, we are free to assign a mass of 1 to the central point representing 
the black hole. Most of the discussion below will be restricted to models with 7 = 1 ("weak 
cusp") and 7 = 2 ("strong cusp"), typical of the values of 7 in bright and faint galaxies 
respectively 

Our assumption of a power-law density dependence with fixed index 7 is reasonable for 
galaxies with strong cusps, in which p ~ r~ 7 , 7 ~ 2, even at radii well outside the sphere 
of influence of the black hole. In weak-cusp galaxies, the shallow inner power law (7 « 1) 
eventually turns over to a steeper dependence at radii r > r b . However, as argued above, 
Tb m r g , and we show below that r g is approximately the radius at which a transition to 
chaos occurs. Thus our 7 = 1 models are expected to yield an accurate description of the 
dynamics of weak-cusp nuclei out to at least the inner edge of the "zone of chaos." 

The gravitational potential corresponding to the stars can be obtained using 
Chandrasekhar's theorem ( Chandrasckhar 1969 , p. 52, theorem 12) which says that for a 



density law that is stratified on similar ellipsoids, the gravitational potential can be written 

as 

*.(*) = -nabcG r W<~)-^" 2 )] dT (3) 

JO ^( r + a 2)( r + 6 2)( r + c 2) 



where 



i"in 2 (T) 

4>{m 2 ) = / p{m' 2 )dm' 2 (4) 

Jrrii. 



2 

ra ax 
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and 



a z + r tr + r c z + r 



For the weak-cusp (7 = 1) case, we have 



<P*(x) = —nabcG / . —dr (6) 

^0 ^(r + a 2 )(T + b 2 )(T + c 2 ) 

while for the strong-cusp (7 = 2) case, 



oc 



^{™Lax) + P° H m max) ~ Po In (A + A + 



$»(x) = -vrafecG / ' " n L^ rfr . (7 ) 



The constant terms, depending on m max , will be ignored in what follows since they have no 
effect on the forces. 

For convenience of numerical calculation, the potential in the strong-cusp case may 
be expressed in terms of a new set of coordinates {r, fi*,u*} which have the following 
definitions (de Zeeuw fc Pfcnniger 1988): 



r 2 = x 2 + y 2 + z 2 , (8a) 

H* = \di + l\[d~ 2 , (8b) 

v* = -di \fd~2, (8c) 

2 2 V K J 



where 



r 2 d x = a 2 (y 2 + z 2 ) + b 2 (z 2 + x 2 )+c 2 (x 2 + y 2 ), (9a) 

r % = [(b 2 -c 2 )x 2 -(c 2 -a 2 )y 2 -(a 2 -b 2 )z 2 ] 2 + 4(a 2 -b 2 )(a 2 -c 2 )y 2 z 2 . (9b) 

In terms of these variables, the stellar potential in the strong-cusp case can be written as 

$,(x)=Alnr + Fi(/i*)+Fi(i/*), (10) 

where 

A = AirGp abcR F (a 2 ,b 2 iC 2 ), (11a) 

r°° lnfr + u) 

FAr) = TiG Po abc / K ' du (lib) 

^0 ^{a 2 + u){b 2 + u){c 2 + u) 
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and 



i2j?(m, n, q) 



2 Jo 



du 



u + m) (u + n) (u + q) 



(12) 



is the Carlson elliptic integral ( Carlson 1988 ). 



Forces may be obtained in analytical form in Cartesian coordinates for both the strong 



and weak cusp cases. By taking partial derivatives of 
are found to be 



where 



2nGabcp 



\/a 2 — 6 2 V a 2 — c' 
2nGabcp Q 



In 



*y 



y/a 2 - bWb 2 - c 2 
2nGabcp 



fx 
A 
-i 



tan _± -i- 



In 



Js 



the weak-cusp force components 



ta11 " (I 



(13a) 
(13b) 
(13c) 



h 
h 

h 
h 
h 
h 
h 
h 



X 



V 'a 2 — b 2 V a 2 — c 2 + abc\ 



b 2 



+ 



-a\x 2 + y 2 + z 2 ) 



= a 2 {{b 2 + c 2 - 2a 2 )x 2 + (c 2 - a 2 )?/ 2 + (b 2 - a 2 )* 2 ) 

+ 2a 2 xV a 2 — b 2 V a 2 — c 2 \J x 2 + y 2 + z 2 

= x 2 (6 2 + c 2 ) + y 2 (a 2 + c 2 ) + z 2 {a 2 + o 2 ) - 2&V + y 2 + z 

= 2yVa? 



b 2 Vb 2 — c 2 J x 2 + y 2 + z 2 
(x 2 b 2 c 2 + a 2 c 2 y 2 + z 2 a 2 b 2 ) - y 2 (b 2 - c 2 )(a 2 - 6 2 ) - b\x 2 + y 2 + z 2 ) 



2y\/a 2 — b 2 \Jb 2 — c 2 J c 2 b 2 x 2 + a 2 c 2 y 2 + a 2 b 2 z 2 



= (zy/ a 2 — c 2 a/6 2 — c 2 + y a 2 c 2 y 2 + b 2 c 2 x 2 + b 2 a 2 z 2 ) 2 — c 4 (x 2 + y 2 + z 2 
= c 2 ((o 2 - c 2 )x 2 + (a 2 - c 2 )?/ 2 + (a 2 + 6 2 - 2c 2 )* 2 ) 
+ 2c 2 zV a 2 — c 2 Vb 2 — c 2 J x 2 + y 2 + z 2 



(14a) 

(14b) 
(14c) 

(14d) 
(14e) 

(14f) 
(14g) 

(14h) 



For the strong-cusp forces, we take partial derivatives of ( PH| ) in Cartesian coordinates 
( |de Zeeuw fc Pfenniger 1988| ) and the force components are given by 



x <9$ a 



r dr 



2x(p-b 2 )(p-c 2 )d$* 2x(u-b 2 )(u-c 2 )d$* 
+ — ~ h — — , (15a) 



u 



dp 



du 
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_ j/gg* | 2 j / (/i-a 2 )(/i-c 2 )^ 2y (z/-a 2 )(z/-c 2 )^ 

r ar r z // — z/ a/i r z u — p, op 

z <9$* 2z (/i - a 2 ) (/i - 6 2 ) 9$, 2^ (i/ - a 2 ) (i/ - b 2 ) <9$* .„ r . 

= o I 2 + — , (15c) 

r or r z fi — v op, r z v — \i ov 



where 



4:irGabcr- l Po R f (a 2 ,b 2 ,c 2 ), (16a) 
nGabcp Rj(a 2 ,b 2 ,c 2 , fi), (16b) 



<9/i 3 
2 

-7rGabcp Rj{a 2 ,b 2 ,c 2 ,i>), (16c 



and 



du 3 

Rj(m,n,q,r) = -/ (17) 

^ (w + rW (w + m)(w + n)(w + g) 



3 A°° 

r) 

is the Carlson elliptic integral 



The magnitude of the radial force in the weak cusp case is constant as a function of 
distance from the center, while in the strong cusp case, the force diverges as r _1 . 

The dynamical time Tu(E) is defined below as the period of a circular orbit of the 
same energy in the equivalent spherical potential, which is defined to have a scale length 
a ave = \J abc. The energy of a circular orbit in the spherical models is 

E c (r) = 3nGp ra ave - (7 = 1), (18a) 

2r 



E c (r) = 2irGp a 
and its period is 

TJr) 



2 

ave 



2 In I — ] -1 



GM, , , , . 

(7 = 2), (18b) 



Gp fa ave y + GM. 



7r(3 — 7) V r J 4ir 2 r 3 



(19) 



Given the energy E in the triaxial potential, we set E c (r c ) = E and solve for r c and T c (r c ). 
The latter is equated to T D (E). 

Equations ( ggg ) and (pb|) are based on the same choices for the zero point of the 
potential as was made in equations @ and ([?]). 

The largest Liapunov exponent was computed for all orbits in the standard way, by 
integrating the equations of motion of an infinitesimal perturbation. Analytical partial 
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derivatives of the forces may be found for both the weak and strong cusp cases. In the 
strong cusp case, the expressions may be simplified using the identity 

dRj(a 2 ,b 2 ,c 2 ,T) 3 1 ( 1 1 _J_ \ p , 2 , 2 2 v 

J i? D (6 2 , c 2 , a 2 ) - — - . /g g (a 2 , c 2 , 6 2 ) - - / . R D (a 2 , b\ c 2 ) (20) 



2(a 2 -r) UK ' ' 7 2(6 2 -r) UK ' ' 7 2(c 2 - r) 
where 

Rn(m,n,r) = Rj(m,n,r,r). (21) 

The equations of motion were integrated using the routine "RADAU" of Hairer & 
Wanner (1996). RADAU is a variable time step, implicit Runge-Kutta scheme which 
automatically switches between orders of 5, 9 and 13. Energy conservation was extremely 
good; energy was typically conserved to a few parts in 10 9 over 100 dynamical times. 



As found also in earlier studies (e.g. [Merritt fc Valluri 1996[ ), a histogram of Liapunov 



exponents of orbits at a given energy evolves toward a characteristic form as the orbital 
integration time is increased. The regular orbits produce a spike at small values of a, 
oTd < 10 -1 ' 1 , whose location moves toward the left roughly as the inverse of the integration 
time. The chaotic orbits produce a broader peak with a well-defined maximum, typically 
at aT D ps 0.2, and a tail that extends almost to the regular orbits. The tail corresponds to 
weakly chaotic orbits that are trapped near regular phase space regions for long periods of 
time. As the integration time is increased, the histogram tends toward two well-separated 
and narrow peaks as the chaotic orbits become increasingly indistinguishable. In what 
follows, the identification of chaotic orbits was based both these histograms and on the 
configuration-space pictures of orbits. 

We henceforth adopt units such that G = a = p = 1. 



3. Orbit Families 

Although the stellar distribution in our models is scale-free, the presence of the black 
hole imposes a scale. We expect the orbital population to change systematically with 
energy, i.e. with distance from the black hole. For each mass model, we defined a grid of 
energy values as follows. We first adopted a set of values M+/M., the ratio of enclosed 
stellar mass to black hole mass : 

logio (^) = {-0.1, 0.0, 0.1, 0.2, • • • 1.7, 1.8}. (22) 
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Each value of M*/M. defines an ellipsoidal surface, m = m*, such that 

M, = 27iabc Po ml (7 = 1), (23) 

M* = 4nabcp m* (7 = 2). (24) 

We then defined the energy E corresponding to this shell as 

$(a;* = am*, 0, 0). (25) 

Table || (weak cusp) and Table || (strong cusp) give x*, E, log (M*/M.), and Tp for the 
mass models with T = 0.5 and c/a = 0.5. 



We followed the standard practice flSchwarzschild 1993| ; [Merritt fc Fridman 199(j| ) 



of defining two sets of initial-condition spaces. Stationary start space consists of initial 
conditions lying on an equipotential surface with zero velocity. X — Z start space consists 
of starting points in the x — z plane with v x = v z = 0. In an integrable triaxial potential, 
stationary start space generates box orbits while X — Z start space generates mostly tube 
orbits. These two start spaces probably contain most of the orbits in reflection-symmetric 
triaxial potentials flSchwarzschild 1993| ). 

As in any non-integrable potential, orbits may be ranked in a hierarchy depending 
on their phase-space dimensionality (Merritt & Valluri 1999). Stochastic orbits fill a 
five-dimensional region and in configuration space populate the entire accessible volume 
within an equipotential surface. Regular, non-resonant orbits occupy 3-tori and densely fill 
some more restricted volume. Resonant orbits satisfy a single relation of the form 

3 

= (26) 

i=i 

between the fundamental frequencies Ui, where the are integers, not all of which are 
zero. Resonant orbits occupy 2-tori and densely fill sheets in configuration space; when 
stable, they have associated with them families of non-thin orbits which mimic the shape of 
the parent thin orbit. Orbits satisfying two independent resonance relations are reduced in 
dimensionality yet again to closed, or periodic, orbits. Periodic orbits are characterized by 
a single base frequency in terms of which the frequency of motion in any coordinate (e.g. x, 
y, z) can be expressed as an integer multiple. 

We label families of orbits associated with a single resonance by the integer vector 
(mi,m 2 ,m 3 ) that defines the resonance. For orbital families associated with doubly 
resonant, or periodic, orbits we use the notation V\ : : z/ 3 , the ratios of the frequencies in 
x, y and z. 
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Figures [I] - [| show the major families of orbits and their starting points in our triaxial 
potentials. Stochastic orbits are present at all energies in both start spaces. They are most 
prevalent in stationary start space, particularly from starting points near the x— and y— 
axes. Motion in the vicinity of the z (short) axis tends to be stable. In X — Z start space, 
stochastic orbits are mostly associated with starting points near the zero-velocity curve, 
or in a few cases with the transition regions between the different families of tube orbits. 
As the energy is increased, stationary start space becomes more and more dominated by 
stochastic orbits; this transition is discussed in more detail in the next section. 

Regular motion in stationary start space is dominated by the pyramid orbits. Sridhar 
& Touma (1999) first described similar orbits, which they called "lenses," in planar, 
harmonic-oscillator potentials containing central point masses. Merritt & Valluri (1999) 
demonstrated the existence of the corresponding 3D orbits in triaxial black-hole potentials. 
Pyramid orbits can be described as Keplerian ellipses with one focus lying near the black 
hole, and which precess in x and y due to torques from the background stellar potential. 
They have a roughly rectangular base whose dimensions are fixed by the amplitudes of 
oscillation in x and y. At low energies, pyramids with a range of shapes exist, having bases 
elongated parallel to both the x— and y— axes. However their major elongation (more 
precisely, the elongation of a symmetrical pair of pyramid orbits oriented above and below 
the x — y plane) is parallel to the z (short) axis. This fact makes pyramid orbits less useful 
than classical box orbits for reinforcing the shape of the figure. 

Close inspection of the pyramid orbits in our numerical integrations reveals many 
resonant pyramid families. Two of the most important are shown in Figure [|: a (3, 0, —4) 
resonance between motion in x and z, and a (0, 6, —5) resonance between motion in y and z. 

As the energy is increased, a 2 : 1 resonance appears in the narrowest pyramid orbits 
lying near the short (z) axis. The opening angle of these "banana" orbits increases rapidly 
with increasing energy, as the stationary point moves along the equipotential surface from 
the short to the long (x) axis. Many of the orbits from the banana family are found to be 
associated with a second resonance. Two such, doubly-resonant familes are illustrated in 
Figure |l|: the (2:3:4) (banana-fish) and (3:4:6) (banana-pretzel) orbits. 

The variation in shape of the bananas with energy is shown in Figure [5]. As discussed 
in the next section, the bananas sometimes persist throughout the "zone of chaos" and 
sometimes disappear, then reappear at high energies, with their major elongation parallel 
to the x-axis. Inside of the chaotic zone, their elongation is counter to that of the figure. 

One resonant family that is apparently not associated with the pyramids is the 
(1, —2, 1) family first discussed by Merritt & Valluri (1999). The major elongation of these 
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orbits is parallel to the intermediate (y) axis; they appear both in the weak- and strong-cusp 
potentials. 

The orbit families in X — Z start space are very similar to those in triaxial potentials 
without black holes ( |de Zeeuw 1985| ; (Schwarzschild 1993 ; [Merritt fc Fridman 1996|) . Tube 



orbits avoid the center due to a primary, 1 : 1 resonance in one of the principal planes 
and are relatively unaffected by the presence of the black hole. The inner long-axis tubes 
are important only in nearly prolate potentials. The most important resonance among the 
tube orbits in our potentials is the 2 : 1 resonance in the meridional plane which generates 
saucer orbits ( [Lees fc Schwarzschild 19"92|) . The saucers appear most prominantly in highly 



flattened potentials. 

We note an interesting feature of the motion in our models. The dynamical roles of the 
long and short axes at low energies are approximately reversed compared to their roles at 
large energies, or in triaxial potentials without black holes. The major families of regular, 
boxlike orbits near the black hole - the pyramids and the bananas - are generated from 
Keplerian ellipses oriented along the short (z) axis, while in triaxial potentials without 
central black holes, it is the long (x) axis orbit that generates the boxes and bananas. 
Similarly, stochastic orbits in our models derive mostly from starting points near the (x — y) 
plane, while in non-singular potentials the instability strip lies near the (y — z) plane 
( poodman fc Schwarzschild 198l|) . This reversal is important because it means that most 
of the regular orbits near the black hole have the wrong elongation for supporting a triaxial 
mass distribution. 



4. Transition to Stochasticity 

The most dramatic effect of the black hole is to induce a sudden change to stochasticity 
in stationary start space as the energy is increased. We investigated this transition as a 
function of 7, c/a and T. Accurate orbital integrations were found to be time-consuming, 
particularly at low energies and in the strong-cusp model. We therefore used the E10K 
supercomputer at Rutgers University to distribute the computations over 64 independent 
processing units. At each energy in each potential, 192 orbits were integrated starting from 
the equipotential surface for a time of lOOTp and the Liapunov exponents were computed. 
Using all 64 processors, the elapsed time for each set of 192 orbits was ~ 10 min for 7=1 
and ~ 30 min for 7 = 2. 

The results are summarized in Figure |^ and [7|. At each energy, the fraction F of the 192 
orbits that were found to be stochastic was computed and plotted. Stochastic orbits were 
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identified both by their location in the histogram of Liapunov exponents at a given energy, 
and by plots of the configuration-space trajectories. While this fraction is not an accurate 
reflection of the fraction of phase-space associated with chaotic motion, the transition from 
regularity to chaos is so sudden that there is no need for a more accurate measure. 

The basic character of these plots is always the same. At low energies, log (M*/M.) < 0, 
the motion is almost completely regular, consisting mostly of pyramid orbits. Starting at 
an energy between log (M*/M») ~ and log (M*/M.) pa 0.5, F increases suddenly to ~ 1 
and remains near unity over a range of energies. Finally, at high energies - log (M*/M.) > 1 
for 7 = 1 and log (M*/M«) > 1.5 for 7 = 2 - F begins to drop and the motion returns to a 
mixture of regular and chaotic orbits. 

The existence of a "zone of chaos" near the centers of triaxial potentials containing 
black holes was first noted by Merritt & Valluri (1999). Based on our more complete set of 
numerical experiments, we can make the following statements about how the properties of 
this zone vary with the parameters of the potential. 

1. For a given triaxiality T, chaos sets in at lower energies in more highly elongated 
models. For instance, for T = 0.5 and 7 = 1, F — 0.8 is reached at log (M*/M«) ~ 0.3 for 
c/a = 0.5 and log (M*/M.) « 0.6 for c/a = 0.8. 

2. The transition from F pa to F ~ 1 takes place more rapidly as a function of 
log (M*/M.) in the more elongated models. 

3. The transition to chaos is interrupted by the appearance of the banana orbits, 
particularly in the more elongated potentials with 7 = 2. For instance, for T = 0.5 and 
c/a = 0.5, the chaotic fraction first increases to F ~ 0.7 at log (M*/M.) = 0.3, then 
decreases again to ~ 0.4 at log (M*/M.) = 0.6 due to the bananas before finally increasing 
to F pa 0.8 at log (M*/M.) > 1. The banana orbits in the most highly flattened models 
(c/a < 0.6) manage to persist throughout the chaotic zone and keep the chaotic orbit 
fraction from reaching 100% 

4. The transition to chaos depends only weakly on triaxiality T for a given elongation 

c/a. 

Figure [7] shows stationary start space as a function of energy for a model 
(7 = 2, c/a = 0.5) in which the banana orbits persist, with an associated family of regular 
orbits, throughout the zone of chaos. In mass models where the bananas disappear, the zone 
of chaos ends with the appearance of a stable resonant family at high energies: typically the 
2 : 3 x — z (fish) resonance for 7 = 1, and the 2 : 1 x — z (banana) resonance for 7 = 2. At 
still higher energies, the orbital populations are similar to those described by other authors 
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flUarpintero fc Aguilar 1998| |Papaphilippou fc Laskar 1998 ; Valluri fc Merritt 1998 ; Waclilin 
|fc Ferraz-Mello 1998|): a complex mix of resonant box orbits, stochastic orbits, and tubes. 



5. Summary 

We have investigated the orbital motion in triaxial nuclei with power-law density 
profiles, p ~ r~ 7 , 7 = (1,2), and central point masses representing supermassive black 
holes. The presence of the central point mass divides the phase space into three radial 
regions. At the lowest energies, the motion is essentially regular. The major orbit families 
are the tubes, the pyramids, and a number of families associated with resonances, most 
prominently the 2 : 1 banana resonance. The pyramid orbits are similar in shape to the 
box orbits of integrable triaxial potentials but have their major elongation parallel to the 
short axis, making them less useful for reconstructing an elongated figure. At intermediate 
energies, the tube orbits persist but the pyramid orbits become increasingly chaotic. The 
transition to a "zone of chaos" occurs rapidly in all of the potentials investigated here, at 
an energy where the enclosed stellar mass is a few times the mass of the central point. 
In the most elongated models with 7 = 2, the bananas can persist throughout the zone 
of chaos; their axis of elongation gradually shifts from the short axis to the long axis of 
the figure. At higher energies, stable resonant boxlike orbits begin to appear in stationary 
start space, generated either from closed orbits like the fish or bananas, or from thin orbits 
corresponding to a 3D resonance. 

Our results are limited in their applicability to the central regions of galaxies where 
the stellar density profile can be approximated as a single power law, p ~ r~ 7 . In bright 
elliptical galaxies and bulges, this is the region within r « r&, the break radius, where the 
central, shallow power law turns over to a steeper dependence at large radii. However we 
argued (Table 1) that is approximately the radius within which the gravitational force 
from the black hole dominates that from the stars; and the results of §4 show that this is 
also approximately the radius of transition to the "zone of chaos" induced by the black 
hole. Thus the onset of chaos in the phase space of real triaxial galaxies should occur at 
approximately the same radius or energy as calculated here. 

Our results highlight two different ways in which central black holes would be expected 
to limit the degree of triaxiality of real galactic nuclei. First, the regular orbits associated 
with motion within the "zone of chaos," the pyramids and the tubes, are mostly poorly 
suited to reinforcing the major elongation of the figure. Second, the black hole induces chaos 
in the motion of filled-center orbits like the pyramids, causing them to occupy a region that 
is rounder than that defined by the equidensity contours of the model. We would therefore 
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expect the degree of triaxiality to be limited inside the zone of chaos by the shapes of the 
regular orbits, and within this region by the lack of regular orbits. These expectations will 
be tested in a future study where self-consistent triaxial models will be constructed. 

This work was supported by NSF grant AST 96-17088 and NASA grant NAG 5-6037. 

REFERENCES 

Bender, R. & Saglia, R. 1999, in Galaxy Dynamics, ASP Conf. Ser. 182, ed. D. Merritt, J. 
A. Sellwood & M. Valluri (ASP: Provo), 113 

Carlson, B. C. 1988, A Table of Elliptic Integrals of the Third Kind. Mathematics of 
Computation, Vol. 51, No. 183, 267 

Carpintero, D. D. & Aguilar, L. A. 1998, MNRAS, 298, 1 

Chandrasekhar, S. 1969, Ellipsoidal Figures of Equilibrium [Dover: New York] 

Crane, P. et al. 1993, AJ, 106, 1371 

Dehnen, W. 1993, MNRAS, 265, 250 

de Zeeuw, P. T. 1985, MNRAS, 216, 273 

de Zeeuw, P. T. & Pfenniger, D. 1988, MNRAS, 235, 949 

Ebisuzaki, T., Makino, J. & Okumura, S. K. 1991, Nature, 354, 212 

Ferrarese, L., van den Bosch, F. C, Ford, H. C, Jaffe, W. & O'Connell, R. W. 1994, AJ, 
108, 1598 

Ford, H. C, Tsvetanov, Z. L, Ferrarese, L. & Jaffe, W. 1998, in IAU Symp. 184, The 

Central Regions of the Galaxy and Galaxies, ed. Y. Sofue (Dordrecht: Kluwer), 377 

Gebhardt, K. et al. 1996, AJ, 112, 105 

Gebhardt, K. et al. 2000, AJ, 119, 1157 

Gerhard, O. E. & Binney, J. 1985, MNRAS, 216, 467 

Goodman, J. & Schwarzschild, M. 1981, ApJ, 245, 1087 



- 15 - 

Hairer, E. & Wanner, G. 1996, Solving Ordinary Differential Equations. II. [Berlin: 
Springer] 

Kormendy, J. 1985, ApJ, 292, L9 

Lees, J. F. k Schwarzschild, M. 1992, ApJ, 384, 491 

Macchetto, F. et al. 1997, ApJ, 489, 579 

Makino, J. 1997, ApJ, 478, 58 

Merritt, D. 1999, in Galaxy Dynamics, ASP Conf. Ser. 182, ed. D. Merritt, J. A. Sellwood 
& M. Valluri (ASP: Provo), 164 

Merritt, D. & Fridman, T. 1995, in Fresh Views of Elliptical Galaxies, ASP Conf. Ser. 86, 
ed. A. Buzzoni, A. Renzini & A. Serrano (ASP: Provo), 13 

Merritt, D. & Fridman, T. 1996, ApJ, 460, 136 

Merritt, D. & Valluri, M. 1996, ApJ, 471, 82 

Merritt, D. & Valluri, M. 1999, AJ, 118, 1177 

Papaphilippou, Y. & Laskar, J. 1998, A& A, 329, 451 

Peebles, P. J. E. 1972, Gen. Rel. Grav., 3, 63 

Quillen, A. C. 1999, in Galaxy Dynamics, ASP Conf. Ser. 182, ed. D. Merritt, J. A. 
Sellwood & M. Valluri (ASP: Provo), 138 

Quinlan, G. & Hernquist, L. 1997, New A, 2, 533 

Quinlan, G., Hernquist, L. & Sigurdsson, S. 1995, ApJ, 440, 554 

Ryden, B. 1999, in Galaxy Dynamics, ASP Conf. Ser. 182, ed. D. Merritt, J. A. Sellwood & 
M. Valluri (ASP: Provo), 142 

Schwarzschild, M. 1993, ApJ, 409, 563 

Sridhar, S. & Touma, J. 1999, MNRAS, 303, 483 

Sambhus, N. & Sridhar, S. 2000, preprint 

van der Marel, R. 1994, MNRAS, 270, 271 

Vallrui, M. & Merritt, D. 1998, ApJ, 506, 686 



-16- 

Wachlin, F. C. & Ferraz-Mello, S. 1998, MNRAS, 298, 22 



This preprint was prepared with the AAS IATeX macros v4.0. 



-17- 




Fig. 1.— 

Major families of orbits in triaxial black-hole nuclei. Each set of three frames shows, from 
left to right, projections onto the (x,y), (y,z) and (x,z) planes, (a) Stochastic orbit, (b) 

Short-axis tube orbit, (c) Saucer orbit, a resonant short-axis tube, (d) Inner long-axis tube 
orbit, (e) Outer long-axis tube orbit, (f) (1, —2, 1) resonant orbit, (g) Pyramid orbit, (h) 
(3,0, —4) resonant pyramid orbit, (i) (0,6, —5) resonant pyramid orbit, (j) Banana orbit. 

(k) 2 : 3 : 4 resonant banana orbit. (1) 3 : 4 : 6 resonant banana orbit, (m) 6:7:8 resonant 

orbit. 
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Fig. 2.- 

Stationary start space for mass model with 7 = 1 (weak cusp) and T = c/a = 0.5. Left 
panels show initial positions on one octant of the equipotential surface; x axis is toward the 
lower left and z axis is up. Frames are labelled by log (M*/M.). Circles represent starting 
points of regular orbits and stars represent stochastic orbits. Colors match the colors of the 
orbit families in Figure 1. Right panels show histograms of Liapunov exponents computed 

over 100 dynamical times. 



-19- 




-0.2 -0.1 0.1 0.2 -2 -1.5 -1 -0.5 

log(crT n ) 




-0.2 0.2 -2 -1.5 -1 -0.5 

log(aT D ) 



Fig. 3.- 

Like Figure 2, for 7 



= 2 (strong cusp). As in the weak-cusp case, pyramid orbits dominate 
at low energies and bananas at high energies. 
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Fig. 4.— 

X — Z start space for mass models with 7=1 (weak cusp), T = 0.5 and three values of 
c/a. Energy is fixed to that of the shell with log (M*/M.) = 0.2. Left panels show initial 
positions in the (x, z) plane. Circles represent starting points of regular orbits and stars 
represent stochastic orbits. Colors match the colors of the orbit families in Figure 1. Open 
circles mark the 1 : 1 closed orbits in the principal planes. Right panels show histograms of 
Liapunov exponents computed over 100 dynamical times. As the elongation of the model is 
increased, stochastic and resonant orbits (e.g. the saucers) become more prominent. 
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Fig. 5 .- 

Banana orbits as a function of energy in five models, each with 7 = 2 and T = 0.5. 

Unstable bananas are not shown. 
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7=1 



7=2 




Fig. 6.— Log(M/Mj Log(M/Mj 

Fraction of chaotic orbits in stationary start space for 7 = 1 (weak cusp) and 7 = 2 (strong 

cusp) models SIS db function of energy. 
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Fig. 7.- 

Stationary start space for the model with 7 = 2, T = c/a = 0.5 as a function of energy, 
through the "zone of chaos." The banana family of orbits persists throughout the chaotic 
zone. The starting points of the resonant banana orbits are shown by the open circles. 



Table 1. Structural Parameters for Two Elliptical Galaxies 



Galaxy 7 r b (pc) M.(M Q ) M*{r b )(M Q ) 

NGC 3379 1.07W 5lW 1.35 x 10 8 ( 2 ) 3.46 x 10 8 
M87 1.26^ 315< 3 > 3.6 x 10 9 ( 4 ) 4.75 x 10 9 



^ebhardt et al. 1996 
2 |Gebhardt et al. 2000 
3 |van der Marel 1994 



1 



Macchetto et al. 1997 
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Table 2. Energy Shells for Weak Cusp Potential, T = c/a = 0.5 
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Table 3. Energy Shells for Strong Cusp Potential, T = c/a = 0.5 
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